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Abstract 

A time- derivative preconditioning algorithm that 
is effective over a wide range of flow conditions from 
inviscid to very diffusive flows and from low speed to su- 
personic flows has been developed. This algorithm uses 
a “viscous” set of primary dependent variables to intro- 
duce well-conditioned eigenvalues and to avoid having 
a nonphysical time reversal for viscous flows. The re- 
sulting algorithm also provides a mechanism for con- 
trolling the inviscid and viscous time step parameters 
to be of order one for very diffusive flows, thereby en- 
suring rapid convergence at very viscous flows as well as 
for inviscid flows. Convergence capabilities are demon- 
strated through computation of a wide variety of prob- 
lems. 

1. Introduction 

Time-marching algorithms are widely used for the 
computation of compressible flows. A major advantage 
of these techniques is that they apply to both invis- 
cid and viscous flows and can be used in conjunction 
with virtually any spatial discretization in all Reynolds 
number regimes. In the past two or three decades, time- 
marching schemes have been widely accepted and ap- 
plied as the method of choice for transonic, supersonic, 
and hypersonic flows. 

In the low subsonic Mach number regime, time- 
marching algorithms do not fare as well. When the flow 
velocity becomes small in comparison with the acoustic 
speed, the time-dependent equations become stiff and 
time-marching methods converge very slowly. These 
difficulties are aggravated if the Reynolds number also 
becomes small. Many problems, however, contain some 
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regions with low Mach numbers while other regions are 
decidedly compressible so that the compressible equa- 
tions have to be used throughout. Consequently, one 
must deal with the stiffness of the equations. Repre- 
sentative problems that contain both compressibility 
effects and low-speed regions can be grouped into two 
classes: high-speed flows with embedded regions of low 
velocity; and low-speed flows with temperature differ- 
ences arising from strong heat addition. Examples of 
such flows are described below. 

High-speed flows with embedded regions of low 
velocity are typified by external, transonic flow with 
embedded low-speed regions near stagnation points or 
by internal flows with low velocity flow upstream of a 
choked area. These embedded regions have little effect 
on convergence when the low-speed region is small, but 
they can dominate convergence when the region is large. 
For example, the low-speed region near the stagnation 
point of an isolated airfoil is seldom large enough to 
affect convergence, but the subsonic flow upstream of a 
strongly converging nozzle may control the convergence 
process. 

Low-speed flows that are compressible because of 
density changes induced by heat addition can be rep- 
resented by problems with surface heat transfer or vol- 
umetric heat addition. The most common example of 
volumetric heat addition occurs in combustion prob- 
lems, but additional problems of interest include ad- 
vanced space propulsion concepts such as laser, solar, 
and microwave thermal propulsion [1] in which elec- 
tromagnetic radiation is used to heat a flowing gas. In 
problems such as these, the equations frequently remain 
stiff over the entire computational domain, so that ef- 
ficiency requirements make it imperative that the stiff- 
ness be dealt with directly. 

Several previous researchers [2-12] have considered 
the low Mach number problem for inviscid flows. The 
present paper is directed towards viscous flows. 

2. Previous Work on Eigenvalue Control 

Two distinct methods have been suggested for con- 
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trolling eigenvalue stiffness at low speeds to enhance 
convergence. The first is to premultiply the time deriva- 
tive by a suitable matrix that scales the eigenvalues of 
the system to the same order of magnitude. The other 
is to use a perturbed form of the equations in which the 
physical acoustic waves are replaced by pseudo-acoustic 
modes. We refer to the former as preconditioning and 
to the latter as a perturbation procedure. The precon- 
ditioning method has the advantage that it provides 
a global solution that is valid at all Mach numbers, 
whereas the perturbation method is valid only locally 
in the regime in which the perturbation is carried out. 

A number of preconditioning studies have been re- 
ported previously. Viviand [2] considered generalized 
preconditioning procedures for a class of hyperbolic sys- 
tems representing the Euler equations and gives specific 
rules that ensure that the preconditioned equations re- 
main well posed. Additional details concerning pre- 
conditioning are presented by Peyret and Viviand [3], 
who are the only ones so far to consider precondition- 
ing in viscous flows. Turkel [4] also discusses precon- 
ditioning with applications to both compressible and 
incompressible flows. Finally, Storti et al. [5] have pre- 
sented some recent theoretical work on eigenvalue con- 
trol that agrees well with Viviand’s findings. References 
2 through 5 are primarily theoretical in nature and con- 
sider the generalized concept of preconditioning in any 
computational regime, not just the low Mach number 
regime with which we are here concerned. Numeri- 
cal studies of preconditioning in the low Mach number 
regime have been reported by Briley et al. [6] for isoen- 
ergetic systems and by the present authors and their co- 
workers [7-9] for the Euler equations. Low Mach num- 
ber convergence enhancement by Briley et al. was lim- 
ited to Mach numbers above 0.05. Our results demon- 
strated that preconditioning provided Mach number in- 
dependent convergence at all Mach numbers between 
0.7 and 10“ 6 , In addition some applications of “time in- 
clining” have been reported by Dannenhoffer and Giles 
[10]. All of these previous works have considered the 
inviscid equations. 

The second general method for eliminating eigen- 
value stiffness is to use a perturbed form of the equa- 
tions of motion. The present authors [8,11] have used 
an expansion of the flow variables in terms of the Mach 
number squared to remove the physical acoustic waves 
and replace them by a set of pseudo-acoustic waves 
whose speeds are comparable to the particle velocity. A 
similar perturbation procedure has also been developed 
by Guerra and Gust afsson [12] based upon expansion 
in Mach number. The perturbation method is effec- 
tive for both viscous and inviscid flows and has been 
implemented for numerous applications [13-15]. Other 
perturbation procedures include the one by Rehm and 


Baum [16], which is specialized toward combustion 
problems and has been applied by several authors in- 
cluding Chenoweth and Paolucci [17] and McMurtry et 
al. [18]. Although these perturbation procedures are 
robust, the nature of the perturbation limits their us- 
age to low subsonic flows. Specifically, the methods 
are not adequate for the transonic flow regime. Our 
focus here is on the preconditioning methods, but we 
take advantage of knowledge gained from perturbation 
methods to develop the preconditioning matrix. 

To date, applications of preconditioning to viscous 
flows have not been reported. We attempted several 
numerical studies for viscous flows using the precon- 
ditioning method of [7,9] and found that it is not ef- 
fective in viscous flows. Accordingly, the purpose of 
the present paper is to develop a method which gives 
reasonable convergence for all Reynolds number ranges 
while keeping its effectiveness in inviscid flows as in the 
previous method. 

3. Problem Formulation 


3.1 Equations of Motion 


The two-dimensional compressible Navier-Stokes 
equations using time-derivative preconditioning can be 
written in the following vector form : 


r dQ.dE,dF_ 
dt d* dy 


= H + L(Q V ) 


( 1 ) 


where T is the preconditioning matrix and will take on 
various forms depending on the preconditioning chosen. 
When T is the identity matrix, we recover the standard 
(nonpreconditioned) equations. The vectors in Eq.(l) 
are 


Q = ( P,pu,pv,e) T 
E = ( pu , pu 2 + p, puv , ( e + p)u) T 
F = (pv,puv,pv 2 +p, (e + p)v) T (2) 

H - (0, 0, —pg, q — pgv) T 
Qv = ( p,u,v,T) T 

In these expressions, the dependent vector Q and the 
inviscid flux vectors E and F take their standard form, 
and H is the source vector, which contains both a vol- 
umetric heat addition rate and a gravitational body 
force. The vector Q v represents the “viscous” variables 
that appear in the diffusion operators. 

The variables in Eq,(2) are defined by standard no- 
tation including density p, velocity components u and 
v, pressure p, temperature T, Mach number M and to- 
tal energy per unit volume e. The coordinate system is 
oriented so that the gravitational body force is in the 
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negative y-direction. Heat is added to the flow by spec- 
ifying the form of the volumetric heat addition function 
q. The gravitational term pgv is the work done by the 
gravity force. The pressure obtained from the perfect 
gas relation is 

P =(7- l )( e -|(« 2 + ^ 2 )) (3) 

where 7 is the ratio of specific heats. 

The differential operator for the viscous terms L is 


low Mach number formulations using perturbation ex- 
pansions for compressible flow [8,11]. The low Mach 
number formulation with this set of variables proves to 
give unconditionally stable results for both inviscid and 
viscous calculations. 

In order to simplify the algebra, first we transform 
the conservative form of Eq.(l) to the non-conservative 
form. We define the non-conservative vector as 

Q = (p,u,v,p) T ( 6 ) 


r± R ± + ±u £ + JL r ± + ±r A. 

dx xx dx dx y dy dy™ 1 ** dx dy y dy 

(4) 


The matrices R XX) R xy , Ry Xy and Ry y are diffusion co- 
efficient matrices that include the viscosity p and the 
thermal conductivity k: 


where the tilde represents the non-conservative vari- 
able. We then premultiply Eq.(l) by the Jacobian P " 1 
= dQ/dQ to obtain 
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In these expressions, Stoke 5 s hypothesis is used for the 
second coefficient of viscosity (A = p). 



3.2 Development of Preconditioning Procedure 


Here the preconditioning matrix is dropped for the time 
being. 

One convenient method for obtaining a new pre- 
conditioning matrix, as well as for transforming the 
(p,u,t>,p) system to the (p, u, v, T) system, is as fol- 
lows. First we subtract the continuity equation from 
the energy equation. This introduces the temperature 
form of the time derivative into the energy equation 
(dp/dt in the continuity equation is replaced by dT/dt 
by means of the perfect gas law.) This temperature 
form of the time derivative also appears well-suited to 
the heat diffusion term. The vector form of this step is 
obtained by premultiplying Eq.(7) by a matrix K \ . 


The preconditioning matrix that was chosen in our 
earlier work [7,9] left the continuity and momentum 
equations in their traditional form but modified the en- 
ergy equation such that time derivatives of />, pu, and 
pv were added. Both stability analyses and numerical 
experiments showed that this preconditioning matrix 
was effective for a wide variety of low Mach number 
inviscid calculations [7,9]. Stability analysis of the full 
Navier-Stokes equations, however, as well as numerical 
experiments shows it is unstable at low Reynolds num- 
bers. Detailed analysis [ 8 ] indicates that this instability 
depends primarily upon the Prandtl number. This de- 
pendence suggests that the modified energy equation 
may no longer be well posed in the diffusion limit. To 
circumvent this difficulty, a new preconditioning matrix 
that is effective both in inviscid and viscous calculations 
is developed in the present paper. 

The basic idea of the new time-derivative precondi- 
tioning procedure is to employ a “viscous” set of depen- 
dent variables ( Q v ) as the primary dependent vector. 
This set of variables is inspired by artificial compress- 
ibility methods for incompressible flow [4,7] as well as 


K 1 ^+K 1 P-\j^ + A!L) = KiP~ l (H+L(Q v )) ( 8 ) 


where the sparse matrix K\ is 


Ki = Diag{l,p y p, 1), 

with nonzero element (4,1) = —yRT (9) 


We then convert the dependent vector Q to Q v by using 
the chain rule, 


KlK2 to +KlP ~ l (lh + T ) = K i p ~ 1 ( H + L (Q'>)) 

( 10 ) 

where is defined as the Jacobian, K 2 = dQ/dQ v . 
Now we precondition Eq.(10) by premultiplying the 
time derivative term by to get 


r -^ i + K ' p ' , (f + %) = k>p-\h + hq,)) 

( 11 ) 
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where the preconditioning matrix T v is defined as 


JM? 

0 

0 

0 ^ 

0 
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0 
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0 J 

(1-7) 
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0 

■ypR/ 


This choice of the preconditioning matrix T v intro- 
duces the artificial compressibility formulation of the 
continuity equation and modifies the energy equation 
to its classical temperature form, while keeping momen- 
tum equations in their standard form. These equations 
are also nearly identical to the perturbation equations 
used extensively [8,11,13-15]. Here, /? M 2 is a scaling 
parameter (M is Mach number), the proper choice of 
which is discussed later. Note that the momentum and 
energy equations now appear to be well-suited for vis- 
cous effects, since they are all of the form 


+ = v 2 * (13) 

where <f> represents u, v, and T, respectively. 

The corresponding time-derivative preconditioned 
system for the conservative form of the equations is 
readily found by premultiplying Eq.(ll) by PK± 1 . 


at 


T- 


. dE dF 
6x + dy 

H + L(Q V ) 

as f=PKi 1 T v , 
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(14) 


(15) 


The eigenvalues of the preconditioned system of 
equations (14) are 


A(f" 1 >l w ) = («,«, 



(16) 


eigenvalues are always positive, while one is negative 
for subsonic flow. 


3.3 Convergence Control at Low Reynolds Numbers 


The stability analysis of Eq.(14) shows that the 
new preconditioning algorithm provides appropriate 
stability for all Reynolds and Peclet numbers. This 
new preconditioning algorithm appears to be at least 
as good as the original preconditioning algorithm for 
inviscid flow and is a dramatic improvement over the 
original one for viscous flows. At low Reynolds num- 
ber flow (below Re=500), however, stability character- 
istics show that amplification factors approach unity, 
suggesting slow convergence rates. The reason for this 
behavior can be understood by considering the viscous 
time step parameter (hereinafter referred to as the von 
Neumann number) <r=pAt/pAx 2 , which becomes im- 
portant at low Reynolds numbers. Control of the CFL 
number alone at low Reynolds number makes the von 
Neumann number so large that the approximate factor- 
ization error in the diffusive terms slows convergence. 
For efficient convergence, we should control both the 
CFL number and the von Neumann number to be or- 
der of unity simultaneously. 

In the present paper, the simultaneous control of 
both the CFL number and the von Neumann number 
is obtained by choosing the scaling parameter k in the 
definition /?. The parameter k is chosen as unity (/?= 
7 RT) to get well-conditioned eigenvalues for the invis- 
cid terms. In low Reynolds number flow, however, we 
use this parameter to specify both CFL and von Neu- 
mann numbers. In the inviscid limit, we choose the 
time step by an appropriate CFL number, 


CFL = 


(u(l + kM 2 ) + c l )At 
2Ax 


(18) 


while in the viscous limit the appropriate viscous time 
step parameter is the von Neumann number. By solv- 
ing Eq.(18) and the above von Neumann number defi- 
nition for k, we find 


where A v =dE/dQ v and the pseudo- acoustic speed c* is 
defined as 


a(a — 1) 
" M 2 (a - 1 + 


(19) 


c ' 2 = u 2 ( 1-^;) 2 + 4PM 2 ( 17 ) 

In order to get well-conditioned eigenvalues for 
the inviscid case, we choose the scaling parameter 
(3=k^RT, where k is a parameter to be determined 
later. For values of k of order unity, this choice of /? al- 
lows the pseudo- acoustic wave speed c f to be the same 
order as the particle velocity u and ensures that three 


where a is CFL f<rRt& x and Re&x is the cell Reynolds 
number (= puAxjp ). Now, by setting the CFL number 
and the von Neumann number to order unity and com- 
puting other variables from given conditions, we can 
obtain the scaling factor k. In general, for a wide range 
;>f Reynolds number, the parameter k can be expressed 


k = Max[ 1, 


a(a — 1 ) 


A/ 2 (o — 1 + 2 ~t-) 


:1 


( 20 ) 


4 



4. Numerical Solution Procedure 

Appropriate preconditioning enhances convergence 
of either explicit or implicit algorithms [9], Here, the 
numerical solution of Eq.(14) is obtained by using an 
Euler implicit discretization in time along with central 
differencing in space. For the efficient solution of the 
resulting matrix, we use an approximate factorization 
such as the Douglas-Gunn procedure [19,20]. This leads 
to 

i 7+A,s ' , <f-£ R -4> ] < 2i > 

[I + AtS ~ 1 (^ - |;«yy|;)]AQ» = -A tS-'R 

where R is the residual of the steady state version of 
Eq.(14) 

flr ftp 

+ < 22 > 

Here A, B and D are Jacobians of the vectors E, F, and 
H and the matrix S is defined as S = f - A ID. This 
formulation differs from the traditional approximately 
factored algorithm only in the calculation of the precon- 
ditioning matrix, and hence additional computational 
cost is negligible. 

5. Boundary Conditions 

In the present study, Method of Characteristics 
based boundary conditions [11] are used at inflow and 
outflow boundaries. In this procedure, it is impera- 
tive to incorporate the preconditioning matrix to re- 
flect the character of pseudo- acoustic waves. To apply 
the Method of Characteristics procedure, we first pre- 
multiply Eq.(21) by the modal matrix containing 
the left eigenvector of the Jacobian r _1 A (or V ~ 1 B) at 
a constant x (or y) boundary. We then multiply the 
result by a selection matrix L that selects those charac- 
teristic equations that represent outgoing information 
at the boundaries to give 

“■'i' + + a - |r»4)]aq, 

= -AtLM~ l S~ 1 Ft (23) 

In Eq.(23), the selection matrix L has different 
forms depending on the boundaries of interest. For in- 
flow boundaries where the flow is subsonic, L becomes 
L = Diag (0, 0, 0, 1), where the nonzero entry repre- 
sents the outgoing characteristic equation (u-c'), while 
the zero entries require physically meaningful boundary 
conditions be specified. For the present study, stagna- 
tion pressure, stagnation temperature, and flow angle 
(v/u) are fixed. 


The outflow boundary conditions are also obtained 
in a similar way. When flow is subsonic, a selection 
matrix is L = Diag (1, 1, 1, 0) and a constant static 
pressure is imposed. For supersonic flow, the selection 
matrix L becomes the identity matrix and no boundary 
condition needs to be specified. 

The boundary conditions imposed on the solid sur- 
face are the no-slip conditions and the normal pres- 
sure gradient condition obtained from the normal mo- 
mentum equation. In addition, either a constant wall 
temperature or an adiabatic wall condition is specified 
depending on the problem. The axis of symmetry is 
treated as a regular field point using symmetry condi- 
tions in lieu of boundary conditions. 

6. Results 

Representative results have been obtained for var- 
ious problems including flow past an isolated airfoil, 
flow through a converging nozzle, and flow in a ther- 
mally driven cavity. In all cases it was demonstrated 
that the preconditioning changed only the rate of con- 
vergence, not the final results. Consequently, we fo- 
cus on comparing the convergence rates of the various 
problems with and without preconditioning. The sav- 
ing realized with preconditioning ranges from a factor 
of two to several orders of magnitude. 

6.1 Flow Past an Isolated Airfoil 

The first test problem considers inviscid and vis- 
cous flow past a NACA0012 airfoil at zero angle of 
attack. A C-type grid (56x31) is used and the outer 
boundary is located 5 chord lengths away from the air- 
foil. Since the flowfield is symmetric, only the half do- 
main is considered. For the wall, a slip boundary con- 
dition is used for inviscid flows, while no-slip, constant 
wall temperature, and zero normal pressure gradient 
are specified for viscous flows. At the outer bound- 
ary, a stagnation pressure, a stagnation temperature, 
and an inflow angle are specified at the inflow region, 
while a constant pressure is specified at the downstream 
end. The remaining conditions come from the Method 
of Characteristics. The outer free stream Mach number 
is 10“ 4 . 

Figure 1 shows the convergence characteristics us- 
ing preconditioning for a wide range of flow conditions 
from inviscid flow to a very viscous flow of Re=l. The 
inviscid case indicates a rate of convergence of one 
order of magnitude per 200 iterations. Convergence 
without preconditioning required some 100,000 itera- 
tions for one order of convergence. For viscous cases 
at Reynolds numbers of 1, 100, and 1000, it can be 
seen that the convergence rate is faster than that of 
the inviscid case by a factor of two. This convergence 
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enhancement may be because physical diffusion serves 
to dissipate the errors in the solution more rapidly. At 
Re= 10000, however, a slight slowdown in convergence 
can be noted. This is because divergence was noted 
when the same CFL number (CFL=6) as that used in 
the low Reynolds number calculations was employed. 
In all cases, it is clear that the preconditioning proce- 
dure enables convergence at these very incompressible 
flow speeds with efficiencies that are competitive with 
those normally observed in subsonic flows. 

The corresponding velocity fields for the three 
Reynolds number cases are shown in Fig. 2. The 

top figures show full domain solutions, while the bot- 
tom figures show magnified results around the airfoil. 
At Reynolds number 1, the flow is nearly symmet- 
ric around the airfoil, approaching Stokes flow. As 
Reynolds number increases, the boundary layer gets 
thinner and at Reynolds number of 10000, separation 
is observed at 80 % chord. The solution accuracy has 
been verified by comparing the pressure coefficient with 
incompressible panel method results for the inviscid 
case but this comparison is not presented for reasons 
of space. 

6.2 Flow Through a Strongly Converging Nozzle 

The second test problem concerns the flow through 
a converging-diverging nozzle. This test problem typi- 
fies a high-speed flow with an embedded region of low 
velocity in which the compressibility effect is signifi- 
cant due to the presence of transonic flow. Two nozzle 
geometries are considered with area ratios (AR) of 10 
and 200. These nozzle geometries are of interest to us 
for applications in solar propulsion [1], which require 
large convergence because of the dilute nature of solar 
radiation. 

First we present convergence characteristics with- 
out preconditioning for AR=1Q for several Mach num- 
bers in Fig. 3. The calculations are made for invis- 
cid flows, and the nozzle geometry is shown in Fig.4. 
Throughout the calculations, an H-grid (71x31) is used 
and the same CFL number is employed. The Mach 
numbers shown in Fig. 3 correspond to the Mach num- 
bers at the throat region. For the M=1.0 case, the 
Mach number variation in the flowfield ranges from 
M=0.05 in the upstream section to M=1.6 at the nozzle 
exit, while for the M=0.1 case, it ranges from M=0.01 
in the upstream section to M=0.08 at the exit. It can be 
noted that from M=0.7 to 1.0 cases, convergence rates 
are nearly the same and it takes 2000 steps to reach 
the 10~ 8 level. This convergence rate is slow at least 
by a factor of seven, compared to the one with precon- 
ditioning shown in the same figure. The convergence 
rate with preconditioning can be understood as the one 
without preconditioning for high subsonic speed flow, 


and thus the convergence comparison suggests that the 
embedded low subsonic flow dominates the convergence 
process and slows the convergence rates. It is also inter- 
esting to note that for the M=1.0 case even though the 
flow passes the sonic line where an eigenvalue stiffness 
exists, the convergence rate is not affected and main- 
tains the same convergence rate as the M=0.9 case. 
This suggests that small regions with the sonic value 
may have little effect on the convergence rates. For 
Mach numbers below 0.5, they gradually start to slow 
as Mach number decreases, as expected. 

Additional calculations for AR=10 and AR=200 
for viscous flows were made. The 71x31 grid for AR=10 
and the 100x50 grid for AR=200 are used, and they are 
highly clustered to the wall. In both cases, calculations 
go through the transonic region. Reynolds numbers 
based on the throat diameter are 6xl0 6 for AR=10 and 
9xl0 5 for AR=200, respectively. Corresponding Mach 
number results using the preconditioning method are 
shown in Figs .4 and 5. Because of the strong nozzle 
convergence, the Mach numbers in the upstream section 
are about 0.05 (for AR=10) and 0.002 (for AR=200), 
respectively. Without preconditioning it is found that 
the solutions are not developed, particularly in the up- 
stream section where the Mach number is low. In Fig. 6, 
the convergence rates with and without preconditioning 
are compared. Without preconditioning, it takes 450 
steps for AR=10 and 2100 steps for AR=200 to reduce 
one order of magnitude in AQ/Q , while with precondi- 
tioning it takes 50 steps for AR=10 and 100 steps for 
AR=200. Thus, the preconditioning method enhances 
convergence rates by a factor of nine for AR=10 and 
twenty one for AR=200, respectively. 

6.3 Flow in a Thermally Driven Cavity 

Our final test problem considers a buoyancy-driven 
flow in a square enclosure. The configuration consists 
of two insulated horizontal walls and two vertical walls 
at temperatures 7* and T c . Even though this prob- 
lem has been a classical natural convection problem 
for many decades, most studies are based on the in- 
compressible formulation with the Boussinesq approx- 
imation [21], which requires a small temperature dif- 
ference between the vertical walls. However, practical 
applications, such as furnace or nuclear reactor design, 
require much larger temperature differences, and com- 
pressible formulations without the Boussinesq assump- 
tion should be employed. Thus, for the present study, 
this problem is chosen as representative of a group of 
low- speed flows that are compressible because of den- 
sity changes induced by surface heat transfer. 

It is known that this problem exhibits com- 
plex flow features depending on the Rayleigh num- 
ber (RaL=p 2 g(3(T h - T c )L 3 C p /fik ), the aspect ratio, 


6 


and the temperature difference parameter (e = (7}» — 
T c )/(Th + T c )). Here /? is the thermal expansion coef- 
ficient and L is the enclosure length. For the present 
study, three Rayleigh number cases, Ra=10 3 ,10 5 , and 
10 6 , are considered with the temperature difference pa- 
rameter € = 0.6. The aspect ratio of this problem is 
one and the transport properties (fi and k ) are evalu- 
ated by using Sutherland’s law. The Prandtl number 
based on reference transport properties is 0.7. A 61x61 
uniform grid is used for Ra^lO 3 , and a 91x91 uniform 
grid is used for Ra=10 5 and 10 6 . 

Convergence rates with preconditioning for the 
three cases are shown in Fig. 7. In all cases, the conver- 
gence behavior shows two distinct slopes. At the initial 
stage, convergence is fast, but it slows down after a 
certain step. The most prominent case is Ra=10 3 . At 
the initial convergence process, 80 steps are required 
for one order of magnitude drop of A Q/Q, but after 
500 steps the convergence rate becomes much slower 
and 1000 steps are required for one order of magnitude 
drop of AQ/Q. The reason for this convergence behav- 
ior is uncertain, but some possible reasons may include 
use of a nonoptimum CFL number or an improper time 
scaling parameter. A further study to understand this 
behavior is necessary. The convergence rate with pre- 
conditioning, however, is fast enough to get solutions 
in reasonable CPU times. 

Figure 8 shows streamline and temperature isolines 
for Ra=10 3 ,10 5 , and 10 6 . It is well known that solu- 
tions with the Boussinesq approximation display a full 
antisymmetric flowfield with respect to the center of 
the enclosure. The present results show a pronounced 
difference and the flowfield is asymmetric. At Ra=10 3 , 
the shift of center of vortex towards the cold wall and 
downwards to the lower wall of the enclosure is man- 
ifest. In all cases, the basic form of the flowfield is a 
recirculating roll. This recirculation is driven by vor- 
ticity generated by the horizontal temperature gradient 
( dT/dx ). At Ra=10 3 , dT/dx is negative over the entire 
flowfield and a single primary clockwise rotating roll is 
formed. At Ra=10 5 and 10 6 , there are two secondary 
rolls embedded in the single roll base flow. These sec- 
ondary rolls are generated because at high Ra numbers 
the intense development of thermal boundary layers in 
the vicinity of the wall leads to the opposite sign of 
dT/dx. Also it can be noted that as Ra increases, the 
secondary rolls intensify and their centers move towards 
the side walls. The accuracy of numerical solutions is 
verified by comparing the Nusselt number at the left 
side wall with a correlation by Chenoweth and Paolucci 
[17] in Fig.9. Good agreement can be observed for all 
three cases. 


Summary 

Extension of the time-derivative preconditioning 
method to viscous flows has been considered. A previ- 
ously developed preconditioning method fails in viscous 
flows because of nonphysical time reversal for diffusive 
terms. In order to circumvent the difficulty, a “viscous” 
set of primary dependent variables is introduced. With 
these primary dependent variables, the equations de- 
generate to the classical diffusion equations in the limit 
of highly diffusive flows. The proper scaling of the time 
derivatives is made to obtain well-conditioned eigenval- 
ues for efficient convergence in inviscid flows. The re- 
sulting algorithm also provides a mechanism for keeping 
both the von Neumann number and the CFL number 
of order one at very viscous conditions, thereby ensur- 
ing rapid convergence at low Reynolds numbers. The 
quantitative effects of the new preconditioning method 
on the convergence of a time-marching algorithm have 
been investigated for several types of problems. It is 
shown that convergence enhancement of the precondi- 
tioning method ranges from a factor of two to several 
orders of magnitude. 
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CONVERGENCE 



Figure 1. — Convergence rates using preconditioning method 
for flow past a NACA0012 airfoil for in viscid and viscous 
cases, M = 10"" 4 . For Re = 1, 100, 1000, CFL = 6, a = 3; 
for Re = 10 000, CFL = 3, k = 1 and for inviscid case 
CFL = 6, k = 1. 
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L2 NORM OF DQ/Q 




Re « 1 

Figure 2. — Velocity contour plots for 


Re = 100 Re = 10000 

viscous flow past a NACA0012 airfoil for Reynolds numbers of 1 , 100, and 10 000. 


CONVERGENCE 



Rgure 3. — Convergence rates with and without preconditioning 
method for inviscid flow in strongly converging-diverging nozzle 
(AR = 10). M t refers to the throat Mach number. CFL = 4, 
k= 1. 



Figure 4. — Mach number contour plot in a strongly converging- 
diverging nozzle (AR = 10). 
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Figure 5. — Mach number contour plot in 
a strongly converging-diverging nozzle 
(AR = 200). 


CONVERGENCE 



Figure 6. — Convergence rates with and without preconditioning 
method for viscous flows in a strongly con verging -diverging 
nozzle (AR = 10 and 200). CFL = 7, k = 1 . 



ITERATION NO. 


Figure 7. — Convergence rates with preconditioning method for 
viscous flow In a thermally driven cavity for Ra = 10 3 , 10 5 
andIO 6 . CFL = 8, k = 10. 


Ra = 10 3 


Ra = 10 s 


Ra = 10“ 


Figure 8. — Streamline and isoline temperature contours for a viscous flow in a thermally driven cavity 


for Ra = 10 3 . 10 5 and 10 6 



Figure 9. — Comparison of Nusselt number with a correlation by 
Chenoweth and Paouilicci [171. 
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